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FOREWORD 

This final report presents work conducted for the Marshall Space Flight Center 
(MSFC), National Aeronautics and Space Administration, in response to the requirements 
of Contract NAS8-38185. The work presented here was performed by REMTECH inc., 
Huntsville, AL, and is titled “External Tank Aerothermal Design Criteria Verification." 

The project manager for this project was Dr. Sarat C. Praharaj. The NASA contract 
monitor for this project was Mr. Lee D. Foster, ED33, of the Induced Environments 
Branch of the Aerophysics Division. 
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ABSTRACT 

If an SSME fails during the initial 1 60 seconds of the Shuttle flight, a return-to-launch- 
site maneuver will be implemented. The period of concern for this task is the pitch-around 
maneuver when the vehicle is flying backward. The intent of this report is to identify and 
define the flowfield at the most critical locations from an environment perspective. The 
solution procedure used to predict the plume heating rates involves both computational 
analysis and engineering modeling. 
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Section 1 
INTRODUCTION 

If an SSME fails during the initial 160 seconds of flight, a return-to-launch-site (RTLS) 
maneuver will be implemented as shown in Fig. 1. The period of concern for this task 
is the pitch-around maneuver where the vehicle is flying backward. Typical altitude and 
angle-of-attack histories during this portion of flight are shown in Figs. 2 and 3, respec- 
tively. Several distinct points have been identified during the turnaround maneuver. 
These points correspond to flight path directions on Fig. 4 where the estimated plume 
boundary and shock location for position F, flying backward, are given. Thus, to define 
the plume with a CFD analysis will require a large grid structure. At positions C and H, 
the free-stream flow along the plume boundary will produce a shear layer which may 
impinge on the ET and orbiter. 

Three flight path positions were chosen for CFD analysis. Currently, the three 
positions which were selected are: 

Position C: High Mach number, possible shear layer impingement along plume bound- 

ary. 

Position F: Flying backward, simplest applied case where plume is axisymmetric, 

similar to validation case. 

Position H: Low altitude, low Mach number, possible shear layer impingement along 

plume boundary. 

The intent is to identify and define the flowfield at the most critical locations from an 
environment perspective. The boundary of the plume will generate a plume shock which 
will affect the airflow over the vehicle, hence, its stability. Plume/air mixture ratios will also 
vary depending upon the location of the vehicle in its trajectory, and this has a strong 
influence on the heat load to the vehicle, specifically to the ET and Orbiter. Heating 
due to aerodynamic frictional heating is dominant if plume impingement is negligible. 
However, a large portion of the hot plume will be forced back to the vehicle due to 
the vehicle’s orientation against the oncoming free stream. Therefore, the most critical 
location from a heating concern is associated with the greatest plume impingement. 
The solution procedure used to predict the plume heating rates involves computational 
analysis (CFD) and engineering modeling (heating indicator approach). 

Following this section, Section 2 details the trajectory used for the RTLS maneuver 
and the modeling (both CFD and engineering) needed to determine the flowfield. Section 
3 discusses the procedure for generating abort environments and provides the environ- 
ments for appropriate body points. Section 4 provided the conclusions reached by the 
Thermal Panel and makes a set of useful recommendations for future work. 
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Return-to-launch-site (RTLS) abort profile is illustrated. Dotted around maneuver is timed for proper main engine cutoff (MECO) 
line shows a normal ascent and normal reentry profile in compari- timing followed by external tank (ET) separation from the orbiter. 

son with the profile that would be flown should an emergency oc- Terminal area energy management phase (TAEM) portion of the 

cur forcing an RTLS during about the first 262 sec of flight. Pitch- abort has tight dynamic loading margins. 

Figure 1 : Return-to-Launch-Site (RTLS) Abort Profile 





Figure 2: RTLS Trajectory: Altitude vs. Time 
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Time Altitude 

A 470 345,254 4.01 

n .. / B 480 347,969 4.00 

Flight Path / c 484 349,198 3.93 



Figure 4: RTLS Flight Path Positions 
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Section 2 

MODELING FOR ABORT CASES 

2.1 Trajectories 

The trajectory used to study the intact-abort return-to-launch-site maneuver was 
provided by Rockwell [1]. For the abort trajectories, the SSME engine #1 out and #3 out 
cases were analyzed. The trajectories were based on statistically acceptable limits which 
would result in conditions to maximize plume heating. Only dispersed trajectories were 
considered in the plume heating calculations. For the above abort cases, la dispersed 
trajectories were used for an SSME #1 out case and an SSME #3 out case. 

Table 1 provides the data for the plume trajectory with SSME #1 engine out (with 
pitot pressure). The power-pitch-around (PPA) portion of the trajectory is marked in this 
table. The trajectory is plotted in Fig. 2. The details of the PPA portion of this trajectory 
are given in Figs. 5 through 7, where altitude, angle-of-attack and Mach number history 
profiles are provided. 


2.2 CFD Modeling 

The following paragraphs describe the CFD study that was conducted to define the 
return-to-launch-site plume-induced thermal environment for the Space Shuttle. The 
task addressed herein is the CFD approach to model the plume flowfield during the PPA 
maneuver. The motivation of the CFD study was to determine the flowfield in sufficient 
detail to perform convective and radiative heating calculations. 

First, it was necessary to develop a CFD procedure by validating the results against 
measured data. The following assumptions were made to simplify the complexities of 
the flowfield: 

1. No yaw was considered. 

2. The flowfields consisted of two frozen species mixing via convection only — one 
frozen specie refers to the exhaust species and the other refers to air. 

3. Diffusion and chemical reactions were neglected. 

The algorithm used for the Navier-Stokes solver consists of the following features: 

1 . Upgrade PARC2D to include the species continuity equation for two species with no 
diffusion and chemical reaction. 

2. Lag the species continuity equation behind the N-S equations. 

3. Couple the species boundary conditions implicitly to species continuity equation. 

4. Transfer information from the species continuity equation to the flow solver via mixture 
gamma and transport properties. 

5. Use the original PARC2D solution as a restart to two-species computation. 
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Figure 5: SSME #1 Out Generic RTLS/PPA Trajectory Altitude Profile 
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Figure 6: SSME #1 Out Generic RTLS/PPA Trajectory Angle-of-Attack Profile 
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Figure 7: SSME #1 Out Generic RTLS/PPA Trajectory Mach Number Profile 
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The second step to the CFD modeling was to use a composite plume for two 
SSME plumes. Then the CFD methodology developed in the first step was used for 
an axisymmetric flowfield corresponding to location F in Fig. 4. This is the trajectory 
point where the angle-of-attack is 1 80 deg. In order to cover the whole range of angle- 
of-attack in the PPA maneuver, trajectory locations varying from C through H (Fig. 4) 
must be covered. At CFD points C and H, for example, a possibility of shear layer 
impingement on the vehicle surface exists. Of these two locations, location C seems to 
be more critical where the angle-of-attack is on the order of -135 deg, and there is a 
possibility of shear layer impingement on the tank base. 

Only the first two steps of the problem were covered during this investigation. The 
two technical notes resulting from the first two steps of this work are documented in 
appendices A and B. 

2.3 Engineering Modeling (Heating Indicator Approach) 

While conducting the CFD study, an easier method has to be sought to estimate 
the upper level of heating on the ET during the power pitch-around maneuver of the 
RTLS. If a conservative method of prediction can be found which produces results not 
impacting the ET TPS design, then flight assessments and design environments can 
be produced in a simple manner for this extremely complex flow process. An early 
attempt was made by Rader [2] in the early 1970’s at Northrop. This method has also 
been utilized by Rockwell in a modified way. Both of these methods disregard most of 
the free-stream/plume interaction effects. The plume gas is assumed to be reversed 
and impinge at free-stream conditions on the ET components producing the worst case 
heating. These heating rates are computed on a stagnation point. 

Review of the method developed by Rader and used in part by Rockwell reveals 
interesting results. The correlation due to Rader is essentially a different version of a 
Fay-Riddell equation, 


q = 



where P * 

Reff 

n 

k 


Stagnation pressure 
Effective radius 

0, 2D 

1, axisymmetric 

Stagnation heating rate plume parameter 


Figure 8 shows the variation of k with p s for various /?«//• For smaller values of R, i.e., 
at R e ff = 0.406 in., the k values are near the frozen flow limit. As R is increased, 
the k values move toward the equilibrium values. These computations were all made 
by Rader [2]. A cross plot of the ratio of k and k eq vs. radius was made in Fig. 8 
for two external values of stagnation pressure. Since no data are available beyond 
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R f f = 8 in., an extrapolation of the plotted data was performed up to R - 10.0 ft. The 
ratio extrapolates to a value of 0.8 or higher, indicating that the equilibrium value is 
being reached for large radii. 

In order to generate a new heating indicator procedure, a conservative philosophy 
was applied at REMTECH. The following assumptions were made in running a heating 
indicator to an RTLS trajectory. 

1. The flow considered chemistry effects. 

Whether the flow chemistry is in equilibrium or the flow has finite reaction 
rates was considered by the code by turning on the appropriate flag in the 
input. In order to check the chemistry of the flow, the equilibrium option 
must be turned on and the results compared against the finite rate results. 

The code used for doing these computations was BLIMPK [3]. 

2. The total enthalpy of the gas, H u was assumed to be equal to total enthalpy of 
the plume gas. 

In other words, no mixture of plume gas and air exists, thus giving the 
most conservative value for H t . 

3. The total pressure of the gas, p t , was assumed to be post-shock total pressure 
based on free-stream air. 

No angle-of-attack effects were considered in the heating indicator cal- 
culations. 
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Section 3 

GENERATION OF ABORT ENVIRONMENTS 

The total scope of this section should cover the environments of various body points 
located on the ET. However, complete environments over the various body points are 
not available now for reasons that will become clear later. The objective of generating 
these environments is to check any exceedances over design. 

Work was initiated to compute the abort environments on the nose of the ET aft 
dome. This point may be the hottest of all other body points located in the rear of the 
ET. The radius of the ET dome (=10 ft) was used to compute a heating indicator with 
the assumptions made in Section 2.3. Two design trajectories, one for aeroheating and 
the other for plume heating, were examined to compute the above heating. The altitude- 
velocity profiles near the PPA maneuver are given in Fig. 9. The heating rate history for 
the aeroheating trajectory is given in Fig. 10. During the PPA period, the highest heating 
occurs around 530 sec with a heating rate of 2 BTU/ft 2 -sec. A similar set of calculations 
was made for the plume heating trajectory. Figure 11 shows the heating rate history. 
Assuming equilibrium edge properties, both equilibrium and finite rate computations were 
made using BLIMPK. The differences between these two computations are minimal. 
Both of these computations used 100 percent plume gas. When compared against 
Rl data book environments, the disparities are obvious. The Rl peak heating is 0.9, 
as opposed to a value of 2.5 from the heating indicator. For comparison, the air-only 
computations are also given in this figure. A number of discussions, teleconferences 
and presentations have been held among REMTECH/MSFC, Lockheed/JSC and Rl to 
resolve these differences. Although a lot of understanding of the methods from different 
organizations has developed, a definite consensus has yet to be reached as to the validity 
of the Rl approach. Figure 1 1 shows results from Lockheed/JSC which are essentially 
the same as these from REMTECH. Another consideration was applied to REMTECH’s 
computations. Since the PPA maneuver occurs in a relatively rarer atmosphere, it is 
worthwhile to check if the heating is transitional. Accordingly, a bridging relationship due 
to Sherman was used in the following way: 

Bridging Relationship: (Sherman) 


N u , 


N, 


UF Af 



where N UFM = free molecular value (obtained from Schaaf and Chambr6 [4] and 
N Uc = continuum value (computed from the BLIMPK code). 

Then ’ ah 

At < = 537.6 sec (plume trajectory), St c = 0.0106 and St pm = 0.089. Stt = 0.00947 
=> q = 2.35 BTU/ft 2 -sec. 


16 


G 





gsg 


2C 




RTR 227-01 



Figure 9: Altitude-Velocity 





Generic RTLS Trajectory (azimuth = 0°) 
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Figure 10: Heating Rate History for Aeroheating Trajectory 
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Figure 1 1 : Heating Rate History for Plume Heating Trajectory 
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The data point was plotted in Fig. 1 1 showing that the transitional heating computation 
did not bring the level of peak heating to the Rl level. 

Since the CFD computations for both 180 deg and 135 deg angle-of-attack have not 
been completed, it was not possible to compute either the convective or the radiative 
environments. Although the 180 deg case has been reported in Appendix B, there seems 
to be some sensitivity to the grid. As a result, this flowfield could not be used for heating 
environments calculations. 
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Section 4 

CONCLUSIONS AND RECOMMENDATIONS 


4.1 Conclusions 

The following conclusions have been made by the thermal panel after the numerous 
discussions and presentations: 

1. The basis for Rl/Downey’s RTLS methodology is not clear. 

2. A heating indicator approach alone will not be able to provide thermal environments 
for design purposes. 

3. However, a heating indicator can provide relative environments from one trajectory 
to the next. 

4. A suitable CFD study is required to define the flowfield for the complex interaction 
region during a PPA maneuver. 

5. The completed CFD studies both by REMTECH and other companies are not accu- 
rate enough to define the required thermal environments over the ET. 

6. Several conclusions were drawn from the CFD experience of Rockwell/Huntsville and 
REMTECH in application to this problem: 

a. The total flowfield is so large that initial gridding is difficult for proper capturing 
of the shock structure. 

b. Gridding for the total flowfield makes gridding for detailed flow around the Shuttle 
configuration difficult. 

c. The subsonic flowfield is quite large, making the placement of the outflow plume 
boundary inaccurate if it is not placed far enough away from the mixing region. 


4.2 Recommendations 

Problems remain in validating the CFD flowfield for the Shuttle PPA maneuver. Since 
no wind tunnel or flight data are available, the CFD codes must be evaluated based upon 
the physics of the flowfield and sensitivity to the grid structure. The required tasks to 

be addressed are: 

Task 1 — Grid Sensitivity Studies 

1.1. Perform a detailed grid sensitivity study of 180 deg angle-of-attack case 

1.2. Exercise the 3-D code for 135 deg angle-of-attack and provide details of 
the plume-air mixtures and thermodynamic properties of the flow near the 
Orbiter/ET base. 
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Task 2 — Physical Phenomena Validation 

2.1. Use available data such as CALSPAN test data as benchmark cases to 
evaluate the modeling of 

a. Hot gas/air turbulent mixing, and 

b. Shock/shear layer interaction 

Task 3 — Document results and prepare a presentation of the results. Give salient 
features, conclusions and recommendations and provide graphical data to 
support each aspect of the problem. 
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Appendix A 

Two-Species Variable Gamma Formulation 
for a Forward-Firing Rocket Nozzle 
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BACKGROUND 

In the event of an SSME failure during the initial 160 seconds of the shuttle flight, a return- 
to-launch-site (RTLS) maneuver will be implemented. The period of concern is the power 
pitch-around maneuver where the vehicle is flying backward at angles-of-attack close to 
180 deg. At these attitudes, a very complex flowfield involving severe free-stream/plume 
interaction affects results. This interaction includes mixing of a high Mach number hot 
exhaust gas from the two SSME’s with the low Mach number free-stream flow. It is 
believed that the hot gas is basically frozen in concentration levels and is allowed to mix 
with air through convection and diffusion. The object of this technical note is to develop 
an algorithm for an axisymmetric case to study the effects of variable 7 on the flowfield 
of a two-species mixture. Since the species diffusion may not play a significant part in 
the solution of the species continuity equation, it has been neglected in this analysis. 


ANALYSIS 


The PARC 2 DR [ 1 ] code was used as the foundation for the Navier-Stokes solver, and 
the species continuity equation for mass fraction, c, of the exhaust species was coupled 
in an explicit fashion to the Navier-Stokes equations via the variable gamma. In two 
dimensions, the species continuity equation, after it was approximated by dropping the 
production term and neglecting the diffusion term, is written in physical space coordinates 
(f, x, y) as 


dc 

dt 


dc dc 
+ u— + v— 
dx dy 


= 0 


( 1 ) 


which is transformed to computational space coordinates (r, f, 77) as 


dc TT dc dc n 
^ + k — + V— = 


/. 


( 2 ) 
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where r -t, £ = £(x, y, t), 7 = 77(1, y, t) and the derivatives are transformed by 


U and V given by 


dt - d T + + TJtdrj 

d x = + rj x d v 

dy = + ly^T] 


U = £ t + 

y = T]t + ur] x + V7] y 


(3) 


(4) 


are the contravariant velocities in the £ and 7 directions respectively, and £ z , £ y , tj x and 
7 y are the metrics etc.^) required for the transformation of the velocities, u and 

v, in the physical space (in the x and y directions, respectively) to computational space. 
The time derivatives of £ and 7 are zero since the grid is stationary. 

Although the Navier-Stokes equations are solved in a steady state mode, the time 
derivative of the dependent variable is retained so that the resulting equations can be 
advanced in time steps of A r to a steady state, thereby allowing the solution to adjust to 
a converged condition over many iterations. Similarly, the species continuity equation (2) 
retains the time derivative and is solved in a steady-state fashion, where each solution 
at time step n is used as the initial condition for the next time step n + 1; this process 
is repeated until the solution c and the solution vector for the Navier-Stokes equations 
[p, pu, pv, pttot ] do not change within a predetermined tolerance. Equation (2) can be 
written in discretized form as 


c n+ 1 - c n + Ar [£/<% + Vd n ]c n+1 = 0 

., . r\ 0 1 

Wlth = 3>r a? 

and ^ = Iri 


• (5) 


Rewriting Eq. (5) in a direction-split form (ADI, or alternating direction-implicit), results in 

[/ + ArUd(:][I + ATVdrj}c n+1 =c n (6) 


in which / is the unit matrix, and a term (A T 2 UVd^d v ), which is a second order term 
in time, has been neglected. At small Ar, this term is much smaller than terms of the 
order Ar. Therefore, Eq. (6) is the required form for the species continuity equation, 
and it is first order accurate in time. 


2 



IS^TTElGhH 


RTN 227-01 


In addition, upwinding is necessary so that communication between neighboring points 
is physically correct. For example, if computations are being performed at grid point 
(j, k), and the 

U — ► 

• • • 

j - 1, k j, k j + 1, k 

contravariant velocity U is in the direction toward (j, k) from (j- 1, k ), then the derivative 
of c in the direction of £ would be 


c £ = c i,k “ c j-hk 


Since cy.jt is unknown and is a known value, upwinding also provides computational 
stability.’ Covering all possible signs for c, upwinding for c in the ^-direction becomes 



and in the j?-direction, 


d n c 


c j,k c j—l,k 

, U> 0 



<3 

II 

o 

(7) 

c j,k ~ c j+l,k 

, U <0 


c j,k ~ c j,k - 1 

, V > 0 


2 

c 

II 

(8) 

°j,k — °j,k+i 

, V'<0 



Equations (7) and (8) are associated with the following form of Eq. (6), where absolute 
values have been taken of U and V: 


[/ + &T\U\d ( } [I + At| V|ajc n+1 = c" (9) 

The solution of Eq. (9) involves solving each direction separately: 

[/ + Ar|V|d I? ]c n+1 = c n (10) 

c" +1 =f (11) 


where 1 

c n = [/ + At|D’|3 { ] V* 

? = [I + &T\V\dJ~ l C" 

= [/+ ArlV’ia,]- 1 [J + Ar|£/ 13 { ] ~ l c" 


( 12 ) 

(13) 
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Inversion of the bracketed terms in Eqs. (12) and (13) is accomplished by first noting that 
the resulting solution at point j, k will involve, at the condition U = 0 (say) points j + 1, k 
and j - 1, k, which corresponds to an equation of no more than 3 points. Therefore, a 
tridiagonal solver is used to invert Eqs. (12) and (13). The value of Ar that is used in 
the species continuity equation is kept equal to the Ar computed from the Navier-Stokes 
flow solver portion of the code. Using the same Ar allows the solution to progress in 
“time” at the same rate in both the flow solver and the mass fraction solver. Furthermore, 
the tridiagonal inversion procedure incorporates the mass fraction boundary conditions 
implicitly to increase the accuracy of and reduce the convergence time of the species 
continuity solver. 

The reason for choosing an implicit, instead of an explicit, mass fraction solver is de- 
scribed briefly below. 

When considering the explicit mass fraction solver, only the interior flowfield points are 
directly involved in the inversion procedure. The boundary points are imposed on the 
interior points explicitly. For example, if a boundary condition is fixed (free stream) or free 
(supersonic outflow), these conditions are set at the points corresponding to the indices 
associated with those boundaries. After this step, the flowfield points are inverted and 
the boundary points are supposed to feed back into the interior through derivatives. 

/ Boundary Conditions 


INTERIOR 

POINTS 

(TRIDIAGONAL 

INVERSION) 


This is the same procedure which is used to invert the flow solver or Navier-Stokes (N- 
S) equations in PARC. However, whereas the interior points in the N-S solver transfer 
information to the boundary points via the derivative term in velocity, there is no such 
transfer from the interior [Eq. (9)] to the boundary for the mass fraction. For example, 
the tridiagonal terms in the matrix below are D, E and F: 
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where E is along the diagonal and corresponds to the points (j,fc )> D corresponds to 
points O' - 1,*) and F corresponds to points (j + 1 ,k). Eq. (9) can be written for the 
; -direction only as 

[1+ArUd^] c n+l = c n (14) 

so that if U is greater than zero, then Eq. 1 4 can be written as 


e^ + ArlT^ 1 -<££*)=«£* 


for a single value of j and Jfe. Equation (15) can be rewritten in matrix form for different 
values of j and constant k as 

r(l + Ar{/) 0 ‘ ' C 2 ,i ' n+1 ' C2,l ’ ” 

(—Art/) (1 + A tU) C3.i c 3.i 


where j goes from 2 to Jm(j = 


(-A tU) (1 + ArCOJ L Om, 1 J 
1 and j = jmax are boundaries), 

Jm = jmax 1 


D = - A tU 
E- 1 + A tU 


F = 0 

As Eq. (15) indicates for j = 2, values of c at location (1, fc) would be required to evaluate 
c 2 However, location (1, ifc) is a boundary point and cannot be included in the vector of 
c values in Eq. (1 6) since this is a matrix inversion for interior points only. Consequently, 
the boundary information for c is not being communicated to the interior points, resulting 
in a non-changing value for C 2 ,*, regardless of the initial values (for U < 0, the same 
argument holds true at inflow boundaries; since there are boundaries in the forward- 
firing nozzle problem in which U < 0 coexists with U > 0, the final value of c will be 
zero everywhere, regardless of the initial values, due to the influence of both types of 
boundaries on the interior point calculations). 

Another way of describing the situation above is by making reference to the N-S equa- 
tions. The explicit solution of the N-S equations as programmed in the original PARC 
code did not encounter such boundary communication problems for the solution to the 
interior points. The corresponding vector c in Eq. (16) would be a vector of Q values 
such that in a single direction, 

[I + A r dfA] A Q n+l = B.HS (19) 
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where A is proportional to U. Since the N-S equations include terms related to the spatial 
derivative of velocity ((<% U) Q n+1 ), the tridiagonal terms, for a boundary point at (j - 1 , k) 
and interior points at (j,Jc) and {j + l,k), are proportional to 

D cc At U(j — 1, k) 

E a 1 -f At C/(j, k) (20) 

F oc At U(j + l,k) 

The situation for the mass fraction equation is such that terms (U d ( c n+1 ) result, leading 

D,F cc At U(j,k ) / 2n 

E oc 1 + AtU (j, k) [ ) 

where the U is now evaluated at point O’, k) only, so that there is no communication 
from the boundary to the interior through U. Although derivative terms for c do appear 
in Eq. (9), these are part of the solution vector c n+1 whose indices only include the 
interior points. 

Therefore, the implicit method was applied to the solution of the mass fraction equation 
by appending the boundary condition equation to the interior point Eq. (9). The form for 
the boundary condition equations is different from those of the interior according to the 
boundary type — fixed (Dirichlet-type), isothermal (Neumann-type), etc. 


Once the mass fraction c of the exhaust species is computed at all points, the ratio of 
specific heats 7 for the mixture is found as 

_ CCy, + (l-e)c K _ Cf^ 

ITTl 1/1 \ 

c + (1 — C)Cy 2 C Vm 

Cp and Cy are the specific heats at constant pressure and constant volume, respectively, 
for species 1 or 2, and subscript m denotes a mixture. Coupling between the species 
continuity equation and the flow solver is accomplished through an equation of state of 
the form 

P = (7m — 1 )P € int (23) 

where p is pressure, p is density and e int is internal energy, with 7 m following from 
Eq. (22). 


To accommodate the difference in the free-stream and exit nozzle value of gamma in 
the normalization of PARC, the following changes had to be made. 


Within the PARC code, the equation of 
free-stream conditions as 

P 

p = 2" 

Poo a 00 


variables are normalized with respect to 


P 


- _ a 

a oo 


state 
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where the overbar denotes a non-dimensional quantity, a is the speed of sound, T is the 
temperature, and «, corresponds to free stream. 

In the calculation of the speed of sound or temperature, the updated code utilizes the 
nondimensional equations: 


2 7 mP 

a = — 

_2 ImP 
— ► a — 

(25) 

p 

P 


p = pRT — 

?P loop 

P 

(26) 


where the 7m in Eq. (25) was equal to 7oo in the original code [1], since 7 was the 
same at every point in the flow. For the present variable gamma approach, the constant 
GAMAREF was changed to GAMMA(J,K) (as computed from Eq. 22), with the reference 
or constant (free-stream) value GAMAFS used in the calculation of temperature only and 
the local value of GAMMA(J.K) used to determine the speed of sound. 

One final consideration is needed to correctly account for the two-species variable 
gamma effects on the flowfield. This becomes important whenever the two species 
have different atomic weights. The equation of state is written in the usual form for a 
single species as 

p = pBT (27) 


where R = and R is the universal gas constant. Eq. (27) can be normalized as 

p PoqUoo P PooR T Too 


or, 

P = P T 

or, 

•e 1 
II 

iT4, 


a 


00 


since 


7co R Too — rigo 

For two species, R must be replaced by the gas constant for the mixture, Rm'. 

p = p R m T where, 

R 


Rm — 


M, 


m 


and M m is the local mixture mass, so that upon normalizing Eq. (30), 

_ —rp Rm Too 

P — pT 


a 


00 


(28) 


(29) 


(30) 


(31) 


7 
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The mixture gas constant involves a summation over each species i, 


B-m — 'y ^ B.\ 

i 

P = pT^Y, c < Ri 

U 00 j 

= V Cj from Eq. (29) 

Too #oo Y 


Also, 



— C Vi 


Substituting Eq. (34) into Eq. (33) results in 


_ 7 ff J2 c i( c Pi ~ C Vi) 

__ pT % 



B. 


00 


(32) 


(33) 


(34) 


= EL SPu — £2sl from Eq. (22) 

Too Boo 


(35) 


= pT_ (t m - l)Cy m 
Too B. oo 

Finally, the equation of state for two species can retain the form it had for a single 
species by writing Eq. (35) as 


P = 



(36) 


where 


T = 


Too goo 
(Tm “ l) c vm 


A similar normalization of the first law of thermodynamics results in 


P = (Tm - 1) P tint 


(37) 

(38) 


in nondimensional form. The temperature is then obtained from Eq. (36) in place of 
Eq. (26). 
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RESULTS 

In order to validate the above algorithm, a few cases were examined. Results for a 
forward-firing sonic jet are shown for two cases below. Both cases have a ratio of 
specific heats equal to 1.382 in the external region of the flowfield and a free-stream 
Mach number of 6.7. Case 1 uses a 7 of 1.4 at the nozzle exit, and Case 2 uses a 7 
of 1.3 at the nozzle exit. Case 1 corresponds to the experimental test case [2], which 
has nitrogen as the exhaust gas and a combusted mixture of methane and air as the 
free-stream gas. In both cases the conditions in the free stream are 

Moo = 6.7 
Poo = 0.3084 psia 
Too = 332.53° R 
T*oo = 3318°R 
Re x /it = 1-4 x 106 
Too = 1.382 

and the total conditions in the sonic nozzle are 

T t = 495° R 

p t - 44.206 psia 
Ptj/Pt 2 = 2 - 46 

Since the O/F ratio of the nozzle exhaust was not known from any of the literature 
available, evaluation of the exhaust and free-stream conditions was made in the following 
manner. Determination of a mixture 7 at each flowfield point requires knowing cp, and 
c*,. Values of cp for the free stream and exhaust species were obtained from JANNAF 
tables, where air was assumed to be the free-stream fluid, and N 2 was taken as the 
exhaust specie. The corresponding value of 7 for air is 1 .382, which was obtained from 
[3], and the exhaust 7 was chosen to be either 1 .4 (Case 1 ) or 1 .3 (Case 2). The specific 
heat at constant volume was then calculated from the relation 


This procedure allows the computation of to be made, but since the value of Rm is 
sensitive to the physically correct values of Cp ( and c^, 7 will not be physically accurate 
(correct values of cp, c„ and 7 can be determined for the exhaust and free-stream gas 
when the ODE (One-Dimensional-Equilibrium) computer code is run). Therefore, within 
the species solver, 7 was taken to be equal to 7 m . Diffusion was assumed to be negligible 
for the blunt-body flowfield of the forward-firing nozzle. Inclusion of diffusion will be 
exercised at a later date to assess the importance of this term. A figure taken from Ref. 
[2] is shown in Fig. 1 to clarify the nomenclature used in the following discussion. 


9 


ren/itech 


RTN 227-01 


Case 1 

Contour plots of pressure, mass fraction of the exhaust species (c = 1 at the nozzle 
exhaust and c - 0 in the external region) and of the ratio of specific heats is shown in 
Figs. 3 through 5, where a plot of normalized pressure contours is given in Fig. 2 for the 
case where 7= 1.4 everywhere for a single species model. 

The value of gamma is equal to 1 .4 everywhere in Fig. 2; in Fig. 3, gamma in the free 
stream is 1 .382 and at the nozzle exit is 1 .4. Comparison of Figs.2 and 3 indicate 
equivalent pressure contour values and shock locations. In addition, the mainstream 
bow shock is slightly closer to the nozzle for the case where 7 in the free stream is 1 .382 
(Fiq. 3) as compared to the case where 7 is 1.4 in the free stream (Fig. 2). Figures 4 
and 5 show similar trends in c and 7, respectively, in which the values at the nozzle exit 
change smoothly to the free-stream values as one moves from the nozzle outward. The 
location of maximum change in c and 7 occurs in a region normal to the nozzle axis and 
is associated with the shear layer located between the mainstream bow shock and the jet 
normal shock. Figure 6 is a side-by-side comparison of Mach number contours generated 
by the variable-species option and the experimental results of [2]. The reattachment 
shock starting behind the jet lip is clearly seen in the Mach number plot. 

The next set of plots depicts the case in which 7* = 1.4, but includes a Wilke mixture 
formulation for both viscosity p and thermal conductivity k to account for the effects 
of a frozen mixture of two species. This is accomplished by utilizing a Sutherland 
formula for each species i to compute individual transport properties m and These 
values are combined along with the mass fraction a into a mixture value for p and k. 
Figures 7 through 10 show no difference in the behavior of the flowfield as compared to 
the Sutherland approximation for a single species as a function of temperature (Figs. 3 
through 6). Nevertheless, the Wilke formulation will be important for problems in which 
diffusion is not negligible and/or problems with more than two dominant species. 


Case 2 

In anticipation of a forward-firing SSME engine with an equivalent mass flow rate equal to 
two SSME’s, whose nozzle exhaust gamma is 1 .3, the gamma for the nozzle in this study 
was changed to 1 .3 (while retaining a free-stream gamma of 1 .382). Results for this case 
appears in Figs. 11 to 14. Pressure contours (Fig. 11) indicate that the jet normal shock 
is somewhat farther from the nozzle exit as compared to the case where gamma at the 
exit was equal to 1.4. A calculation of p*lp t for a sonic nozzle (sonic divided by total) 
shows that the pressure is indeed larger at the nozzle exit when the value of gamma is 
smaller The shock system between Case 1 and Case 2 is very similar. Exhaust mass 
fraction contours (Fig. 12) range from 1.0 at the nozzle exit to 0.0 in the free stream 
with a similar trend as in Case 1. The region of lower contour values near the nozzle 
corner is due to recirculation which can be seen in the Mach number plot. This region of 
reduced exhaust mass fraction occurs in the expansion fan emanating from the jet lip. 
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Evidently, lower value of c corresponding to the regions farther from the nozzle made 
their way to the expansion fan and did not reduce to zero as it did in Case 1 . 

Figure 13 indicates the gamma contours with 7 = 1.3 at the nozzle exit and 7 = 1.382 in 
the free stream. The region of largest gradient is along the axis and corresponds to the 
jet-free-stream interface or shear layer. As the flow travels to the backside of the nozzle, 
the (convective) mixing region extends over a larger portion of the flowfield; the gradients 
in the mass fraction and gamma are more gradual on account of the decreased velocity. 
A Mach contour plot (Fig. 1 4) describes the velocity field as in Case 1 , where the jet 
normal shock stands farther out from the nozzle exit with the lower 7 . The secondary 
reattachment shock is clearly indicated in Fig. 14. 

The final plot (Fig. 15) is a side-by-side comparison of the Mach number contours from 
the two-species variable gamma cases (with -y e = 1 .4 and 7 * = 1 .3). This figure reveals 
the difference in the locations of the respective jet normal shocks. Only the Sutherland 
transport property formulation for a single species is used. The lower gamma at the 
exhaust results in a jet normal shock that lies slightly farther from the nozzle exit. 
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Figure 2: 










6: Mach Contour Comparison Between Two-Species Variable Gamma Results (Exhaust 
Gamma = 1 .4, Free-Stream Gamma = 1 .382) and Ref. 2 
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Abstract 

A procedure for solving an axisymmetric SSME nozzle exhausting with an equivalent 
mass flow of two nozzles is presented for the return-to-launch-site (RTLS) maneuver 
for an engine-out case. The algorithm incorporates a two-species variable gamma 
formulation along with diffusion terms associated with a species conservation equation 
Coupling of this algorithm is accomplished explicitly to the Navier-Stokes PARC2DR and 
PARC3DR codes. Axisymmetric results indicate a very complicated flowfield with a bow 
shock extending to approximately 700 radii in front of the nozzle, and low temperatures 
on what would be an external tank (ET) aft of the nozzle. 


Background 

In the event of an SSME failure during the initial 160 seconds of the shuttle flight, a 
retum-to-launch-site (RTLS) maneuver will be implemented. The period of concern is the 
Dowered-pitch-around (PPA) maneuver where the vehicle is flying backward at angles- 
of-attack close to 180 deg. At these attitudes, a very complex flowfield involving severe 
free-stream/plume interaction effects results. This interaction includes mixing of a high 
Mach number hot exhaust gas from the two SSME’s with the low Mach number free- 
stream flow. It is believed that the hot gas is basically frozen in concentration levels and 
is allowed to mix with air through convection and diffusion. The object of this technical 
note is to develop an algorithm for an axisymmetric case to study the effects of vanable 
7 on the flowfield of a two-species mixture. 


Analysis 

Species Solution Procedure 

The PARC 2D R [1] and PARC3DR codes were used as the foundation for the Navier- 
Stokes solver, and the mass fraction, c, obtained from the species continuity equation 
(with and without diffusion) was coupled in an explicit fashion to the Navier-Stokes 
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equations via the variable gamma and diffusion terms if diffusion is considered. In two or 
three dimensions, the species continuity equation, after it was approximated by dropping 
the production term, is written in physical space coordinates (f, x, y) as 

^ + V • (piV) = V • {pDn V c.) W 

at 

where Du is the binary diffusion coefficient, /?,• and c,- are the density and mass fraction of 
species i respectively, and V is the velocity vector. Rewriting Eq. (1) using the definition 
of mass fraction 


p 


results in 


+ V • {p Vci) = V • {pD 12 V ci) ( 3 ) 

at 

wjth v _ d_ j + d_ j + Eq. (3) can be reduced to the following form with the use 
of continuity 


P 


Dei 

~Dt 


= V • (P D n 


Vc.) 


(4) 


where the global continuity equation was used to simplify Eq. (4) and ^ is the substantial 
derivative. 

At this point Eq. (4) can be rewritten in a strong conservation form for three dimen- 
sions as 

a (pc,) d(pua) gfcvci) + g(£wcO = v . {pDuVCi) (5) 

dt dx dy dz 


where u, u and w are the velocities in the x, y and z directions, respectively. Either 
Eq. (4) or Eq. (5) can be used to compute the mass fractions. The reason for using the 
weak form instead of the strong conservation form will be discussed next. 

Upon normalizing the variables according to reference values 


P U - V - P - 1 J7- J. j — t re t 

-i— u= , V = ,p= ' x ~ 1 ’P — / ,' 1 ~ l l t 

p ref a re f a re f IrefPref ‘ re/ l rcf l ref 


( 6 ) 


where a is the speed of sound, Eqs. (4) and (5) can be written as 



dtitCj) 

ct 


+ 4(P « + J:(/> u «) + $-(p w a) - 




(7) 
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Pref a ref [r e/ 


C * 

Sc = ■■ ■— 

P D 12 

The Schmidt number, 5c, will be taken as a variable quantity in this development. 

To transform Eqs. (6) and (7) from physical space to computational space (£, 77, C. 
r) start with the functional relationship between the coordinates: r -t, i - i ( x,y,z,t ), 
t; = 77 (x,y,z,t), < = C so that the derivatives from physical to computational 

space become 

dt = d T + it m d v + C< a c 

= 6 + 7* ft? + Ci ft; (8) 

ft, = £y ft + 7y ft? + Cy ft 

ft = 6 ft + 7z ft + C* ft 

and the contravariant velocities are defined by 

U = it + «£* + 

V = T)t + UT] X + VTJy + WT]z (9) 

w = C< + u Cz + u Cy + tl>Cz 


where ft = |i , etc., are the metrics required for the transformation of the velocities u, v, w 
in the physical space to the computational space. The time derivatives of i, 77 and < are 
zero since the grid is stationary. 

The weak form of the species continuity equation in the computational space be- 


comes 


Cr + U C{ + V Cq + W Cf = 



( 10 ) 


where the overbars and subscript, t, have been dropped for convenience, and the 
derivatives are given by Eq. (8). In this form, the species conservation equation can 
be discretized and the left-hand side can be rewritten in direction-split or alternating- 
direction-implicit (ADI) form and the right-hand side is kept explicit: 


[I + Arf/ft] [I + ArFft] [i + ArWft]c n+1 ’ = RHSD (11) 


where 


RHSD = c" + Ar 


1 

R& r ef 



(12) 


I is the unit matrix, 9 S = etc., and terms of second order in time have been 
neglected. At small A r, these second-order terms are much smaller than terms of order 
Ar. Equation (11) would be advanced in time using current values at time step n as 
initial conditions for the next time stop n+ 1 until steady-state convergence was achieved. 
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Similarly, the strong form of the species conservation equation can be transformed to 
the computational space after multiplying Eq. (7) by the inverse Jacobian J 1 , resulting 
in the delta form 

[/ + A tU ( ] [I + A tV„\ [I + A r W c ] A (^) " +1 = RHS + RHSD (1 3) 



An additional term arising from this process 
can be shown to be equal to zero [2]. 

The decision as to which of the diffusion equations to be solved — the weakly 
conserved form of Eq. (1 1 ) or the strongly conserved form of Eq. (1 3) — was made based 
upon the type of numerical stability scheme desired. Solving the strong conservation form 
of the species equation would require either central differencing that initially includes 
artificial dissipation (since Eq. (13) is the same form as the fluid dynamic, central 
difference equations in PARC), or upwinding. It was decided that since the species 
continuity equation can be quite easily upwinded, then Eq. (11) would be discretized 
accordingly. 

One-sided differencing is utilized to bias the discretizing species continuity equation 
in the direction of the fluid flow; upwinding provides correct communication between 
neighboring points. For example, if computations are being performed at grid point (j, k), 
and the 

U -* 

• • • 

j-l,k j,k j + l,k 

contravariant velocity U is in the direction toward O’, *0 fro™ O’ - 1> k )< then the derivative 
of c in the direction of £ would be 


c £ = c i,k ~ c j-l,k 


Since c jk is unknown and Cj _ hk is a known value, upwinding also provides computational 
stability as it is naturally dissipative. Covering all possible signs for c, in two dimensions, 
upwinding for c in the f— direction becomes 


Cj, k - Cj- hk 

, tf>0 


Cj+\,k -Cj-XJt 
2 

, U = 0 

(15) 

c j,k ~ Cj+ l.Jfc 

, U< 0 
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and in the ^-direction, 


' c i,k ~ 

, V>0 

2 

< 

II 

o 

V Cj.Jk “ c J,k + 1 

, V<0 


(16) 


Equations (1 5) and (16) are associated with the following form of Eq. (1 1 ), where absolute 
values have been taken of U and V : 


[I + Ar|tf|^][J + Ar|V|d„]c n+1 = RHSD n 


(17) 


The solution of Eq. (1 7) involves solving each direction separately: 

[I + Ar\V\d v }c n+l =7ffiSD n 

n 

c n+1 = RHJD 


(18) 

(19) 


where 


RHSD" = [/ + Ar|tf|0 c ] V 

RHSB = [I + RHSD 

= [/ + ArlVI^f 1 [I + &T\U\ds\~ l c n 


( 20 ) 


( 21 ) 


Inversion of the bracketed terms in Eqs. (20) and (21) is accomplished by first noting that 
the resulting solution at point k will involve, at the condition U = 0 (say) points j + 1, * 
and 7 — 1 k t which corresponds to an equation of no more than 3 points. Therefore, a 
tridiagonal solver is used to invert Eqs. (20) and (21). The value of A r that is used in 
the species continuity equation is kept equal to the Ar computed from the Navier-Stokes 
flow solver portion of the code. Using the same Ar allows the solution to progress in 
“time" at the same rate in both the flow solver and the mass fraction solver. Furthermore, 
the tridiagonal inversion procedure incorporates the mass fraction boundary conditions 
implicitly to increase the accuracy of and reduce the convergence time of the species 
continuity solver. The reason for choosing an implicit, instead of an explicit, mass fraction 
solver is described briefly below. 

When considering the explicit mass fraction solver, only the interior flowfield points 
are directly involved in the inversion procedure. The boundary points are imposed on 
the interior points explicitly. For example, if a boundary condition is fixed (free stream) 
or free (supersonic outflow), these conditions are set at the points corresponding to the 
indices associated with those boundaries. After this step, the flowfield points are inverted 
and the boundary points are supposed to feed back information into the intenor throug 

derivatives. 
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,/ Boundary Conditions 


INTERIOR 
:: POINTS Ilil 
Ilf (T R I D I AGO NAlill' 
INVERSION) 


This is the same procedure which is used to invert the flow solver or Navier-Stokes 
(N-S) equations in PARC. However, whereas the interior points in the N-S solver transfer 
information to the boundary points via the derivative term in velocity, there is no such 
transfer from the interior [Eq. (17)] to the boundary for the mass fraction. For example, 
the tridiagonal terms in the matrix below are D, E and F: 



matrix for interior points in j-direction 


where E is along the diagonal and corresponds to the points (j, k), D corresponds to 
points (j - 1, k) and F corresponds to points (j + l.Jfc). Equation (17) can be written for 
the j-direction only as 


[/ + Ar Uds] c n+1 = RHSD n 


(22) 


so that if U is greater than zero, then Eq. (22) can be written as 



- <£;,*) - r hsd 


},k 


(23) 
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for a single value of j and k. Equation (23) can be rewritten in matrix form for different 
values of j and constant k as 



where j goes from 2 to Jm (j = 1 and j = j ma x are boundaries), 



J m — jmax 1 

(25) 

and 

D = —A tU 



E = 1 + Art/ 

(26) 


F = 0 



As Eq. (23) indicates for j = 2, values of c at location (1, k ) would be required to evaluate 
c 2 fc. However, location (1, k ) is a boundary point and cannot be included in the vector of 
c values in Eq. (24) since this is a matrix inversion for interior points only. Consequently, 
the boundary information for c is not being communicated to the interior points, resulting 
in a non-changing value for C 2 ,jt» regardless of the initial values (for U < 0, the same 
argument holds true at inflow boundaries; since there are boundaries in the forward- 
firing nozzle problem in which U < 0 coexists with U > 0, the final value of c will be 
zero everywhere, regardless of the initial values, due to the influence of both types of 
boundaries on the interior point calculations). 

Another way of describing the situation above is by making reference to the N-S 
equations. The explicit solution of the N-S equations as programmed in the original 
PARC code did not encounter such boundary communication problems for the solution 
to the interior points. The corresponding vector c in Eq. (24) would be a vector of Q 
values such that in a single direction, 

[/ + Ar d^A] A g n+1 = RHS (27) 

where A is proportional to U. Since the N-S equations include terms related to the spatial 
derivative of velocity ((d s U) Q n+l ), the tridiagonal terms, for a boundary point at {j - 1, k) 
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and interior points at (j, k) and (j + l,k), ar© proportional to 

D <x At U(j — 1, fc) 

E oc 1 + At U ( j , k ) (28) 

F oc At U(j + 1, k ) 

The situation for the mass fraction equation is such that terms (U c n+1 ) result, leading 

t0 D,F oc At U(j,k) (29) 

E oc 1 + A TU{j,k) 

where the U is now evaluated at point (j, k) only, so that there is no communication 
from the boundary to the interior through U. Although derivative terms for c do appear 
in Eq. (17), these are part of the solution vector c n+1 whose indices only include the 

interior points. 

Therefore, the implicit method was applied to the solution of the mass fraction 
equation by appending the boundary condition equation to the interior point Eq. (1 7). The 
form for the boundary condition equations is different from those of the interior according 
to the boundary type — fixed (Dirichlet-type), isothermal (Neumann-type), etc. Along a 
solid boundary the mass fraction gradient is set to zero, while along a free boundary the 
mass fraction is either extrapolated from the interior flow as in the subsonic case or set 
equal to the value from the preceding grid point as in the supersonic case. 

Energy Equation 

The contribution of diffusion to the energy equation appears in the heat flux term 
(neglecting the Dufour Effect) 

Vi (30) 

i 

where hi is the enthalpy of component i, and the mass diffusion velocity Vi is related 
to Fick's Law of Diffusion [3] as 


Normalizing Eq. (30) results in the following term in physical space 

1 1 \ 

where the Schmidt number, Sc, is defined by 

c a 
Sc = — t— 

pDij 

The Schmidt number is allowed to vary in this formulation. 


(31) 


(32) 


( 33 ) 
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The diffusion coefficient was computed from the following equation obtained from the 
bifurcation approximation of Bird [4] 


Da = 


D 

FiFj 


(34) 


where D = 1 - 719x ^-— (r) 1 - 659 cm 2 / sec with p in atmospheres and T in degrees Kelvin, 

and Fi = (mif 26 J) 0 - 489 , where m< is i th the species molecular weight. In discretizing 
D\ 2 , it is advisable to smooth out the significant nonlinear variations in its value across 
adjacent cells. This reasoning, plus the differences in the cell sizes, suggests an average 
diffusion coefficient [5] which guarantees the flux entering an interface from the left of a 
discretized cell is equal to the flux exiting that interface at the right. 

n £j £*±iixi±i -X») (35) 

2 A+l (x i+ i - X.j + Di [ X i - Xi+I ) 

where the species subscript 1 2 has been removed as it is understood that D is the binary 
diffusion coefficient, and where I + 1 refers to the actual computational grid center and 
i + i indicates a cell interface: 


i,j + 1 




» + 


Additional Considerations 

Discretization of the species second derivative terms in computational space is 
accomplished according to whether the derivative is mixed or not mixed. For similar 
second-order derivatives, second-order central spatial difference operators are used, 
whereas, when applying a mixed derivative to c, average derivative values are computed 
at locations between the two grid points in a particular derivative direction. 


= (C;+l,Jb,i — c j,k,l) — ( c j,k,l c j-l,k,l) 

c J+ i t jb+i,r ~ Cj+i.i-i,/ , Cj,k+U ~ c j,*-U 
2 2 

1 [ c i,ife+l,r “ c j,k-l,l . c j-l,k+l,l ~ c ;-l,t+l,r 

2 [ 2 2 



( 36 ) 


9 



REIVITECH 


RTN 227-02 


Since the implicit equations are solved or inverted for the a along a single j or k or / line 
while changing one of the other indices (as opposed to solving for the species for the 
entire computational space in a single inversion), provisions were made to incorporate 
cyclic iteration between the possible sweep direction ordering. 

Finally, to provide for axisymmetric diffusion, the species continuity 


dci dci dci 

p ^ + pu d^ + pV ^ 



(37) 


in which v is the radial velocity and y is the radial coordinate, was transformed to 
computational space coordinates and coded in the upgraded PARC code, hereafter 
referred to as PARCSPEC. 

Specfes/Flow Solver Coupling 

Coupling of the results from the species concentration solver to the fluid dynamic 
solver is accomplished in an explicit fashion. The continuity and momentum equations 
remained unchanged in form; however, due to the use of Wilkes mixture formula to 
calculate the transport properties, the energy equation changes its form to incorporate 
a variable thermal conductivity and Cp. The mass fractions, along with the species gas 
constants, are used to compute 7 at every point in the flowfield, which replaces the 
constant 7 required in the original version of PARC. The ratio of specific heats 7 for the 
mixture is found as 

_ c Cpi + (1 — c)cp, _ Cp^ (38) 

^ C C Vi + (1 c )c»2 

C« and c„ are the specific heats at constant pressure and constant volume, respectively, 
for species 1 or 2, and subscript m denotes a mixture. Coupling between the species 
continuity equation and the flow solver is accomplished through an equation of state of 
the form 


P = (7m - l)pCint 


(39) 


where p is pressure, P is density and e int is internal energy, with 7m following from 
Eq. (38). 

To accommodate the difference in the free-stream and exit nozzle value of gamma 
in the normalization of PARC, the following changes had to be made. 

Within the PARC code, the equation of state variables are normalized with respect 
to free-stream conditions as indicated in Eq. (6) plus the additional terms 


T = 


T 

Tr*f 



(40) 


where T is the temperature, and a is the speed of sound. 

In the calculation of the speed of sound or temperature, the updated code utilizes 
the nondimensional equations: 


a 


2 


ImP 

P 



(41) 
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p = pRT — + T = ( 42 ) 

where the 7 m in Eq. (41) was equal to 7 «, in the original code [1], since 7 was the 
same at every point in the flow. For the present variable gamma approach, the constant 
GAMAREF was changed to GAMMA(J.K) (as computed from Eq. 38), with the reference 
or constant (free-stream) value GAMAFS used in the calculation of temperature only and 
the local value of GAMMA(J.K) used to determine the speed of sound. 

One final consideration is needed to correctly account for the two-species variable 
gamma effects on the flowfield. This becomes important whenever the two species have 
different molecular weights. The equation of state is written in the usual form for a single 

species as 

p = pRT ( 4 3) 


where R = £, and R is the universal gas constant. Equation (43) can be normalized as 


or, 


or, 


P Prtf a ref = P Pref R T T rt} 
- -7f TrcfR 

p = pT 


a rtf 


_ PT 

P = z — 

7r cf 


(44) 


since 


7 ref R Tref = “rtf 


(45) 


For two species, R must be replaced by the gas constant for the mixture, Rm. 

P = P RmT where, 

R 


Rm — 


m m 


and m m is the local mixture mass, so that upon normalizing Eq. (46), 

_ _7t7 Rm Tref 

t=pT T 7 ' 

The mixture gas constant involves a summation over each species i, 

Rm ~ 'y ^ Ci r* 


(46) 


(47) 


(48) 
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00 i 

= ^ V Ci Ri from Eq. (45) 

Inf Rref Y 


Also, 


Ri — Cpi — Cv; 


Substituting Eq. (50) into Eq. (49) results in 

-7p 'll Ct'C 0 ?. _ *”»•) 

- = pJL j 

" Iref Rref 


P T C Pm ~ C V. 

7re/ Rref 


from Eq. (38) 


(49) 


(50) 


(51) 


_ P 7* (7m 
7re/ Rref 

Finally, the equation of state for two species can retain the form it had for a single 
species by writing Eq. (51) as 


P = 



(52) 


where 


A similar normalization 


7oo Roo 

(7 m — 

of the first law of thermodynamics results in 


P = (7m - 1) P e int 


(53) 

(54) 


in nondimensional form. The temperature is then obtained from Eq. (52) in place of 
Eq. (42). 


Results 

In order to validate the above algorithm, a few cases were examined. Results for a 
forward-firing sonic jet [ 6 ] for two cases 7 # = 1 - 4 and 7iv = 1 -3 at a 7 f s = 1 .4 compared 
very well with experimental results. Diffusion was not invoked in the code in that study. 

The engine-out condition on a PPA maneuver in an RTLS trajectory was modeled 
as a single nozzle with a mass flux equal to twice that of one SSME, and with a radius 
scaled to V2 times a single SSME radius. A grid was created with INGRID [7] assuming 
the nozzle to be a simple outflow plane, with a straight afterbody with a radius of 1 .05 r ; 
the location of flow impingement on the external tank (ET) (which is not part of the grid 
model) would then be approximated by examining the solution in the region where the 
ET would exist in relation to the SSME’s. 
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Conditions on the nozzle and free stream are as follows: 

Table 1: Nozzle and Free-stream Conditions 



Nozzle 

Free Stream 

p (Ibm/ft 3 ) 

3.226 x 10~ 3 

9.0328 x 10 -6 

T (°R) 

2193.0 

409.18 

p (lbf/in 2 ) 

2.619 

1.3681 x 10 -3 

m (Ibm) 

14.132 

28.858 

Mach Number M 

-4.67 

1.464 

7 

1.20 

1.40 

The free-stream conditions were taken as the reference conditions in PARCSPEC. 


This study focuses on the possible impact on the thermal environment caused by 
the SSME firing in the same direction as the Shuttle flight direction. At a flight angle a 
of 180 deg, the plume impingement on the ET may not occur; however, at a flight angle 
of 135 deg, the free-stream air may force back the plume which is expanding from the 
nozzle into a shear layer that impinges the ET. 

Axlsvmmetrlc SSME 

The extent of the grid ranges from 1454 (two-nozzle equivalent) radii in front of the 
nozzle and 300 radii behind the nozzle, exit plane along the axis, and 2800 radii in 
the radial direction. A grid size of 140 points (axially) and 340 points (radially) was 
constructed with clustering near the axis and nozzle expansion region (Fig. 1). The 
nozzle is located at the x = 0 location on the figure. The large size of the gnd prevents 
the nozzle from being seen on this scale. Due to the very large physical extent of this 
arid the grid lines are highly compacted in the vicinity of the nozzle, and more sparse as 
the grid extends far away from the body. Smoothness between adjacent gnd cells was 
adhered to as best as possible to limit the nonlinearities associated with computations 
across such cells. Diffusion was neglected in this set of runs. 

Convergence was difficult to obtain when the conditions specified in Table 1 were 
imposed at once. Instead, the nozzle conditions were taken at their specified values, 
except that 7iV was set to 1.4, and the free-stream conditions were adjusted gradually 
between restart files until the specified conditions were reached. Each set of free-stream 
conditions was run to a reasonable convergence before adjusting the conditions for the 

next casG. 

Contour plots of streamlines are shown in Figs. 2-3 for the case where 7 jv = 7 /s 
= 1.4 and the molecular weight m N = m f s = 28.858. The demarcations on the axis 
specify a range of 200 radii. A smooth variation in the streamlines is indicated in these 
figures, where the flow turns back at a point equivalent to 180 radii in front of the nozzle 

on the axis. 
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The next step in the convergence process was to reduce m and ms from 1 .4 and 
28.858, respectively, to 1 .2 and 1 4.1 32, which corresponds to the H 2 - 0 2 system value. 
Results for the PPA conditions at a = 180 deg are presented in Figs. 4-8. The bow 
shock is located at roughly 700 radii in front of the nozzle while the nozzle shock occurs at 
250 radii before the nozzle (Fig. 4). Due to the high nozzle pressure and low atmospheric 
conditions, the nozzle shock extends a great distance away from the nozzle and is highly 
curved as a large expansion from the axis to the free-stream condition must occur. The 
values of gamma (Fig. 5) indicate a smooth variation in the extent of the mixing region. 
Along the axis, supersonic Mach numbers of Fig. 6 show the increase of the value of 
4.87 at the nozzle to 10.9 about 170 radii before the nozzle and then a decrease to 
zero at an X/R value of 200 where the two opposing flows stagnate. The behavior of 
the streamlines (Fig. 7) in front of the nozzle is due to the small gamma and molecular 
weight of the exhaust gas compared to that shown in Fig. 2. Clearly, the hot plume gas 
does not impinge on what would be the ET. A close view of the nozzle exit temperature 
profiles (Fig. 8) indeed verifies the low temperature in the region just behind the exhaust. 
In fact, the temperature is close to the free-stream value. 


Results In Progress 

When the nozzle gamma was changed from 1.4 to 1.2, it was found that the 
streamlines fluctuated over a large distance along the outflow plane. This is reasoned 
to be caused by flowfield and exit plane coupling which could be eliminated if the exit 
plane were extended farther downstream. An extension of 400 radii or a total of 700 
radii — comprise this new grid, and results are being generated at the current time. 

In addition, the diffusion upgrade to the two species methodology is being checked 
out for the two-dimensional sonic jet to ascertain the effect that diffusion has on the 
variable gamma solution. 

To study the three-dimensional effects on the RTLS flowfield, the axisymmetric grid 
was revolved about the axis and points between the outer boundary and the bow shock 
were eliminated, as these points were originally employed to capture the unknown 
location of the shock. Eighteen planes comprise a total space of 1 80 deg with half 
of the domain being a mirror image of the constructed grid domain. Computations are 
in progress for this case. 
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Figure 1 : Grid 
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Figure 4: Pressure Contours for 7^ = 1.20, 7 /s = 1.40 
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Figure 6: Supersonic Mach Contours for 7,v = 120, 7 js = 1.40 
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